Tuumn 


L*J  IL’Jll  I  ■  n  * 


1*.  REPORT  SECURITY  CLASSIFICATION 
UNCLASSIFIED 


2».  SECURITY  CLASSIFICATION  AUTHOR 


3.  DISTRIBUTION /AVA, LABILITY  OF  REPORT 
Approved  for  public  release; 
distribution  unlimited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

AFOSK-TK*  8  9  -  1  10  9 


6b.  OFFICE  SYMBOL  |7a.  NAME  OF  MONITORING  ORGANIZATION 
(If  applicable)  f 


,  ^  ^  c  /  l  \  i  *  l  >  v_  /0  j  •  /  VJ 

"j  Ls>  Cc 


7b.  ADDRESS  (City,  State,  and  ZIP  Code) 
Building  410 

Bolling  AFB,  DC  20332-6448 


3a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 


8c  ADORESS  (City,  State,  and  ZIP  Code) 

Building  410 

Bolling  AFB,  DC  20332-6448 


8b.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable) 

NM  Arose-  8?-c\ 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO. 


(mb aapi 

fT««  iyi  [•]  l  ■  k’  [•! 


1 1 .  TITLE  (Include  Security  Classification) 


i  T  i  |4  A  /)  ■  ' 


'  A  1  .  ,.A,  T  C 


61102F 


\J0  l  Lj  f  jL  L  I  i j  > j. l.  'L  i  pT  A ~i  ,  l'  Jvi  O  t  t  L  w  i C  i i  0  (4  5  cl  0 J  I-)  "I  1?  1 1 

1  h  t  L.  l  i  V  •  1  •  t- J  A  1  A  (7  (-  k .  I  C  , ,  S  ■)  C  0  \  1 


14.  DATE  OF  REPORT  (Year,  Month,  Day)  115  PAGE  COUNT 


16.  SUPPLEMENTARY  NOTATION 


COSATI  CODES 


GROUP  SUB-GROUP 


18.  SUBJECT  TERMS  ( Continue  on  ravers*  if  necessary  and  identify  by  block  number) 


19.  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

This  effort  supported  the  research  activities  of  20  researchers  during  their 
visit  to  ICASE,  as  a  result,  10  papers  have  appeared  on  Issues  related  to 
parallel  computation  Including  such  titles  as  "Reordering  computations  for 
parallel  execution",  "Multiprocessor  L/U  decomposition  with  controlled 
fill-in,  and  "Analysis  of  a  parallelized  nonlinear  elliptic  boundary  value 
problems  solver  with  applications  to  reacting  flows". 


20  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

0 UNCLASSIFIED/UNLIMITED  (fit SAME  AS  RPT  □  OTIC  USERS  UNCLASSIFIED 


1 22  b.  TELEPHONE  (Include  Area  Code)  |  22c.  OFFICE  SYMBOL 


00  Form  1473,  JUN  86  Previous  editions  are  obsolete. 


Q  O. 


ICASE  REPORT  NO.  87-63 

AhOtftK-  nt.  8  9-1109 

ICASE 

POLYNOMIAL  APPROXIMATION  OF  FUNCTIONS  OF  MATRICES  AND  ITS 
APPLICATION  TO  THE  SOLUTION  OF  A  GENERAL  SYSTEM  OF  LINEAR 
EQUATIONS 


Hillel  Tal-Ezer 


Contract  No.  NAS1-18107 
September  1987 


INSTITUTE  FOE  COMPUTER  APPLICATIONS  IN  SCIENCE  AND  ENGINEERING 
NASA  Langley  Research  Center,  Hampton,  Virginia  23665 

Operated  by  the  Universities  Space  Research  Association 


NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

I  lampion.  Virginia  23665 


POLYNOMIAL  APPROXIMATION  OF  FUNCTIONS  OF  MATRICES 
AND  ITS  APPLICATION  TO  THE  SOLUTION  OF 
A  GENERAL  SYSTEM  OF  LINEAR  EQUATIONS 


Hillel  Tal-Ezer* 

Raymond  and  Beverly  Sackler  Faculty  of  Exact  Sciences 
Tel- Aviv  University 
Israel 

and 

Institute  for  Computer  Applications  in  Science  and  Engineering 


Abstract 

Frequently,  during  the  process  of  solving  a  mathematical  model  numerically,  we  end 
up  with  a  need  to  operate  on  a  vector  v  by  an  operator  which  can  be  expressed  as  /(A) 
while  A  is  N  x  N  matrix  (ex:  exp(A),  sin(A),  A-1)*  Except  for  very  simple  matrices,  it  is 
impractical  to  construct  the  matrix  /(A)  explicitly.  Usually  an  approximation  to  it  is  used. 
In  the  present  research,  we  develop  an  algorithm  which  uses  a  polynomial  approximation 
to  /(A).  It  is  reduced  to  a  problem  of  approximating  f(z)  by  a  polynomial  in  z  while 
z  belongs  to  the  domain  D  in  the  complex  plane  which  includes  all  the  eigenvalues  of 
A.  This  problem  of  approximation  is  approached  by  interpolating  the  function  /  (z)  in  a 
certain  set  of  points  which  is  known  to  have  some  maximal  properties.  The  approximation 
thus  achieved  is  “almost  best."  Implementing  the  algorithm  to  some  practical  problems  is 
described.  fl~.  - 

Since  a  solution  to  a  linear  system  Ax  =  b  is  x  =  A-1  b,  an  iterative  solution  to  it  can 
be  regarded  as  a  polynomial  approximation  to  /(A)  =  A-1.  Implementing  our  algorithm 
in  this  case  is  described  too. 
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I.  Introduction. 

Let  A  be  a  N  x  N  matrix  whose  elements  belong  to  C,  and  f(z)  is  a  function  such 

that 

/(*):  C-C  .  (1.1) 

The  matrix  f(A)  can  be  defined  in  the  following  way:  If 


w 

f(z)  =  Y2a  i(a-*o)*  \z-zo\<r 


then 


t=0 


oo 


f(A )  =  a*(-A  “  *o/)*  iA*  ~  *>l  <  r 


(1.2) 


(1.3) 


«=o 


where  A*  is  an  eigenvalue  of  ^4.  In  (1.3),  f[A)  is  expressed  as  an  infinite  polynomial.  It  can 
be  shown  [Gant  59]  that  f(A)  as  defined  above  can  be  written  also  as  a  finite  polynomial 
of  degree  <  iV  —  1  as  follows: 

f(A)  =  ^2  [f(Xk)Hki(A)  +  •  •  •  +  /^-x)(Afc)^ (A)]  (1.4) 

where 

A*  -  an  eigenvalue  of  A 

jk  -  multiplicity  of  A*  (1.5) 

Hkj (A)  —  polynomial  in  A  of  degree  <  N  —  1  . 

In  many  practical  applications,  we  would  like  to  compute  a  vector  u>  which  results 
from  operating  with  the  matrix  f(A)  on  a  vector  v 

u  =  [f(A))v  .  (1.6) 

Using  expression  (1.4)  for  this  purpose  has  two  major  disadvantages: 

(a)  The  exact  knowledge  of  the  eigenvalues  is  req  '  d. 

(b)  /(A)  is  expressed  as  a  polynomial  of  degree  <  N—  1  which  can  be  large;  thus  it  results 
in  an  highly  time  consuming  operator. 

In  this  paper,  we  would  like  to  present  an  algorithm  which  approximates  the  matrix  f(A) 
by  a  polynomial  Pm(A),  where  m  <<  N.  Since  (1.4),  we  have 

Sm(A)  =  f(A)  -  Pm(A)  =  £  f£(Ak)Hkl(M)  +  •  ••  +  E^-'\\t)Htii[A)]  (1.7) 

I _ I 


fc»l 


_  Codes 

(Avail  and/or 
Dlst  |  Spoclal 


r 


□□ 
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while 

E(z)  =  /(,)  -  P„(.)  .  (1.8) 

Hence,  our  problem  is  reduced  to  a  problem  of  approximating  the  function  f(z)  by  a 
polynomial  Pm(z)  while  z  belongs  to  a  domain  in  C  which  includes  all  the  eigenvalues  of 
A.  It  is  obvious  from  (1.7)  that  in  the  case  where  jk  >  1,  Pmk~^  has  to  approximate 
fUh-i)  at  the  point  A*. 

In  Section  2,  we  describe  briefly  how  to  approach  this  problem  of  approximation.  Using 
the  results  of  this  section  we  can  construct  a  computational  algorithm,  once  a  suitable 
conformal  mapping  function  is  known.  In  very  few  cues,  it  is  possible  to  get  this  function 
analytically.  Usually  we  have  to  resort  to  a  numerical  method.  A  lot  has  been  done  in  this 
area  of  conformal  mapping  and  there  are  suitable  routines.  In  our  work,  we  have  used  a 
set  of  routines  written  by  N.  Trefethen  based  on  his  paper  [TrefSO],  which  are  very  efficient 
and  reliable.  The  interested  user  can  get  the  routines  from  the  author  upon  request.  The 
algorithm  which  results  from  making  use  of  the  approximating  polynomial  is  presented  in 
Section  3.  In  Section  4,  we  deal  with  the  rate  of  convergence  of  the  suggested  method.  Few 
examples  of  mathematical  models  for  which  the  new  algorithm  can  be  implemented  are 
given  in  Section  5.  The  numerical  solution  of  a  system  of  linear  equations  Ax  =  b  where 
A  is  a  general  nonsymmetric  matrix  can  be  regarded  as  a  particular  case  of  our  problem 
where  the  function  which  has  to  be  approximated  is  f(z)  =  1  jz.  Section  6  is  devoted  to  this 
subject.  Using  tools  from  the  theory  of  approximation  in  the  complex  plane  for  inverting 
nonsymmetric  matrices  has  been  studied  already  by  other  researchers.  A  few  of  them 
are:  T.  Manteuffel  [Mant77],  [Mant78],  D.  Smolarski  and  P.  Saylor  [SmSa85],  Gutknecht 
[Gutk86],  Y.  Saad  [Saad87],and  others.  The  algorithm  discussed  in  our  paper  gains  its 
simplicity  from  the  fact  that  it  is  based  on  the  powerful  tool  of  interpolation.  Thus,  it  can 
be  implemented  in  a  straightforward  way,  once  the  conformal  mapping  function  is  known. 
An  important  conclusion  of  the  analysis  is  that  for  the  general  nonsymmetric  case,  the 
significant  factor  which  governs  the  rate  of  convergence  of  most  of  the  iterative  algorithms 
is  not  the  well  known  condition  number  |)A||  •  || A-1 1|.  It  is  shown  that  the  relevant  factor 
is  p/R  where  p  is  a  parameter  which  measures  the  size  of  the  domain  and  R  is  related  to 
the  orientation  of  the  domain  with  regard  to  the  singular  point  of  the  function  j  (z  =  0). 
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Based  on  this  conclusion  it  is  shown  that  some  iterative  algorithms,  like  standard  SOR, 
are  based  on  an  incomplete  optimization  process.  We  conclude  this  paper  in  Section  7  by 
giving  some  numerical  examples. 


2.  Polynomial  Approximation  in  the  Complex  Plane 

Let  D  be  a  bounded  closed  continuum  in  C  such  that  the  complement  of  D  is  simply 
connected  in  the  extended  plane  and  contains  the  point  at  oo.  Define  now 

A(D)  -  the  space  of  functions  which  are  analytic  at  every  point  of  D. 
nm  -  space  of  complex  polynomials  of  degree  <  m. 

Then  it  is  well  known  [SmLe68]  that  for  every  /  €  A(D)  there  exists  6  jrm  such 

that 

||/-PA!)oo<||/-Pm||oo  VPmG7Tm.  (2.1) 

From  algorithmic  point  of  view,  it  is  relatively  complicated  to  find  P^.  In  many  cases,  it 
is  quite  simple  to  find  a  polynomial  approximation  which  is  “almost”  as  good  as  P^.  It  is 
found  by  projecting  A(D)  on  irm.  If  Sm  is  a  projection  operator 

Sm  :  A[D)  ,  (2.2) 


then 


/  -  Sm(/)  =  f -Pm  +  Smi.Pm  ~  />! 


(2.3) 


thus 

!l/-5m(/)||<(l  +  ||5m||)|!/-P^||.  (2.4) 


If  ||Sm||  is  reasonably  small,  regarding  5m(/)  as  “almost”  as  good  as  P„  is  justified.  For 
example,  if  D  is  the  unit  disc  and 


Sm{!)  = 


/W(0) 


t 


(2.5) 


||Sm|l  ^  ~2  log(m)  +  1  +  0(1)  . 


then  (GeMa75| 


(2.6) 
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If  D  =  [—1,  lj  and 

m 

Sm{f)  =  YlakTk(z)  (2.7) 

k=0 

while 


Tk(z)  —  coa(k(coB~1  z))  (Chebyshev  polynomials)  ,  (2.8) 

then  the  bound  on  ||Sm||  is  exactly  the  same  as  in  the  previous  case  (2.6).  A  generalization 
to  an  arbitrary  domain  D  (as  defined  in  the  beginning  of  this  section)  is  given  by  a  Faber 
polymomials  expansion.  Let  4>[z)  be  a  conformal  mapping  from  z  to  u  which  maps  the 
complement  of  D  to  the  complement  of  a  disc  of  radius  p  such  that 


lim  M 

Z—*OQ  Z 


=  1  . 


(2.9) 


p  is  the  logarithmic  capacity  or  the  transfinite  diameter  of  D  [SmLe68]. 
If  the  Laurent  expansion  at  oo  of  f<£(z)jm  is 


(^(z)]m  ~  Z™  +  cm-lZm  1  +  1  •  •  CiZ  +  Cq  +  C_xZ  1  -f  •  *  •  ,  (2.10) 

then  Faber  polynomial  of  degree  m,  Fm(z),  which  corresponds  to  the  domain  D  is  the 
polynomial  part  of  (2.10): 


We  have 


’m  +  1  + - +  Co  ■ 

(2.11) 

1  f 

*/+* 

(2-12) 

while  R  is  choeen  sufficiently  large  so  that  D  is  contained  in  the  interior  of  the  region 
bounded  by  the  circle  |z|  =  R  [SmLe68].  Given  /  €  A(D)  then 


/(*)  ~  Y  *kFk{*)  • 

k=0 


(2.13) 


.,  =  J -  / 

2*»  JM=r  uk+1 


The  coefficients  a*  are 
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where  ip(u)  is  the  inverse  of  4>(z)  and  r  >  p  is  sufficiently  small  that  /  can  be  extended 
analytically  to  Ir.  (Ir  is  the  closed  Jordan  region  with  boundary  Tr  where  Tr  is  the  image 
under  rj)  of  the  circle  |u/|  =  r.)  In  (Ella83j  it  was  shown  how  the  c}' s  and  a*’s  can  be 
computed  efficiently  using  F.F.T.  Based  on  (2.13),  the  matrix  f(A)  can  be  approximated 
by  Pm(A)  where 

m 

Pm(A)  =  £a*F*(A).  (2.14) 

fc= o 

When  D  is  an  ellipse  with  axes  parallel  to  the  real  and  imaginary  axes, then  Fk  is  the 
translated  and  scaled  Chebyshev  polynomial  Tk  [Mark77j;  thus  it  satisfies  a  simple  three 
term  recurrecnce  relation.  This  recurrence  relation  enables  us  to  compute 


ui  = 


I 


U=0 


(2.15) 


efficiently.  A  detailed  description  of  these  cues,  approximating  the  operator  eA,  can 
be  found  in  [Tale85],  (Tale86j,  [KoTa86],  [TaKo84|.  For  more  complicated  domains,  F* 
satisfies  a  k  terms  recurrence  relation  [Ella83];  thus  from  memory  point  of  view  it  is 
impractical  to  use  an  algorithm  based  on  (2.15). 

This  major  drawback  can  be  overcome  by  using  a  different  approach  to  the  approxi¬ 
mating  problem.  It  uses  interpolation  as  the  projection  operator  Sm.  Using  this  tactic  we 
are  confronted  with  the  following  question:  Which  are  the  interpolating  points  in  D,  such 
that  the  interpolating  polynomial  will  be  “almost”  as  good  as  P^l  It  is  known  [SmLe68] 
that  if  D  is  the  unit  disc,  then  two  possible  sets  of  points  are: 

(a)  Zj  —  0  j  as  1, . . . ,  m  (zeroes  of  zm) 

(b)  The  m  roots  of  unity. 

In  a  similar  way,  for  a  general  domain  £?,  a  two  possible  sets  of  points  are  [SmLe68j: 


(a)  the  m  zeroes  of  Fm{z)  (2.14a) 

(b)  Zj  =  0(wy)j1’  =  1, . . .  ,m,  while  u>j  are  the  m  roots  (2.14b) 

of  the  equation  wm  =  p. 

In  the  first  case,  it  can  be  shown  [GeMa75j  that 

^  ^r!og(m)  + 1  +  0(1) 


(2.15) 


t' 


6 

while  in  the  second  case 

<  - log(m)  ■+•  1  +  0(1)  .  (2.16) 

7T 

Since,  interpolating  at  Zj  =  0(u/j)  j  =  0, . . . , m  does  not  involve  computing  Faber  poly¬ 
nomials,  it  is  the  simplest  and  most  efficient  approach  to  this  problem  of  approximation. 
The  interpolating  polynomial  in  Newton  form  is: 

Pm(z)  =  oo  +  ax(z  -  zo)  +  o2(z  -  20)(z  -  zi)  -(-  •  •  •  +  am(z  -  20)  •  •  •  (z  -  zm_-)  (2.17) 

where  a*  is  the  divided  difference  of  order  k 

f\z o,  •  •  * ,  sjt]  k  =  0, . . . ,  TYi  .  (2.18) 

When  p  ( the  logorithmic  capacity  of  D)  is  large,  o,  will  be  very  small;  thus  it  is  preferable 
to  make  the  following  change  of  variables 

z  —  zj p  .  (2.19) 

Hence 

f{z)  =  f{z)  =  f(p  ■  z)  (2.20) 

and 

Pm(z)  3  /(£)  (2.21) 

where 

Pm(z)  =  bo  +  bi(z  -  z0)  +  •  •  ■  +  bm(z  -  z0)  ■  ■  ■  (z  -  zm-i)  (2.22 a) 

bk  =  f[zo,---,Zk\  k  =  0, . . . , m  .  (2.22 b) 

The  only  difficulty  in  finding  Pm{z)  is  to  get  the  function  t/>(w).  For  simple  domains, 
it  is  possible  to  find  this  function  analytically.  For  more  complicated  domains  one  has  to 
resort  to  a  numerical  approach. 

When  the  domain  D  is  a  polygon,  the  mapping  function  is  Schwartz-Cristoffel  trans¬ 
formation.  In  [TrefSO],  a  very  reliable  and  efficient  algorithm  for  mapping  the  interior  of 
the  unit  disc  to  the  interior  of  the  polygon  is  described.  Since,  in  our  case,  the  mapping 
of  the  exteriors  is  needed,  the  routines  should  be  modified  accordingly. 
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Separated  domains:  In  some  very  important  physical  phenomenon,  we  have  a  situation 
where  the  eigne  values  of  the  matrix  A  are  clustered  in  two  or  more  separated  domains  which 
can  be  far  apart  (for  example:  stiff  O.D.E.’s  problems  of  parazitic  waves  in  meteorological 
sciences,  spectral  methods  for  nonperiodic  boundary  values  problems,  etc.).  Hence,  the 
domain  D  is  a  union  of  k  subdomains: 

D  =  u kimXDi  . 

In  this  case,  the  complement  of  D  is  not  simply  connected  any  more  but  just  connected. 
The  basic  theory  regarding  the  simply  connected  case  extends  to  the  more  general  one.  A 
detailed  analysis  of  this  problem  is  carried  out  in  our  next  paper.  It  is  shown  there  that 
the  interpolating  points  can  be  taken  as  a  union  of  sets  of  points  achieved  by  considering 
each  domain  separately.  In  Section  7,  we  bring  some  numerical  examples  of  this  case. 

3.  The  Algorithms 

Based  on  (2.22),  we  will  approximate  the  operator  f[A)  by  Pm{A)  while 
Pm{A)  —bo I  +  bi(A  -  zol)  +  6a(A  -  z0/)(A  -  h /)+ 

(3-1) 

•  ••  +  6m(o  -  zol)  ■■■{A  -  Zm- ll)  ■ 

A  =  (1/p)  A 

Operating  with  Pm(A)  on  a  vector  v  is  carried  out  by  the  following  algorithm 
Algorithm  1: 

u  =  v 
u  —  bov 

For  t  =  1, . . . ,  m  do  (3.2) 

u  =  [A  -  z,-i/)u 
w  =  w  4-  biU 

The  output  of  Algorithm  1  is  the  vector  w  which  satisfies 

w  =  Pm{A)v  .  (3.3) 

Roundoff  errors  of  Algorithm  1  depend  strongly  on  the  arrangement  of  the  interpolating 
points. 
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In  Appendix  A  we  describe  how  to  arrange  the  points,  taking  into  account  this  phe¬ 
nomenon  of  nrvdoff  errors. 

In  man/  practical  problems,  we  have  the  following  situation: 

1)  A  is  a  real  matrix  (3.4a) 

2)  /(*)  =  7f*I-  (3.4b) 

In  this  case,  it  is  possible  to  design  an  algorithm  similar  to  (3.2)  which  will  be  carried 
out  in  real  arithmetic  even  so  Z{  and  bi  are  complex  numbers.  This  result  is  based  on  the 
following  two  theorems: 

Theorem  X:  Let  zo>Zi, z3,I2,. •  •  ,  Zfc.Zfc  be  2k  interpolating  points  where  zq  and  z\  are 
real  numbers.  If  Pu,-i{z), 


p2k-i{z)  =  Oo  +  fliz  H - a2k- \Z2k  1 


(3.5) 


is  the  interpolating  polynomial  of  a  function  f(z)  which  satisfies  (S. 4b)  then 
ao.ai, . . .  ,a2fc_x  are  real  numbers. 

Proof.  The  Lagrange  formula  of  Pik-x{z)  is 

*  k 

=  E  ijWM  +  E  CW/ft) 

;= o  j~2 

where 


(3.6) 


Qiki*})  z~zi 


j  =  0, . . . ,  k 

(3.7) 

3  =  2,...,* 

(3.8) 

(3.9) 

1=0  is] 


Since  (3.9)  we  have 


k 

n 

t=2 


q%h (z;)  =  n (*j - *<) n (* j  - *tj  j=o,...,k 


(3.10) 
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Qlkfe})  —  II  (7y  2»)  fr  (*/  -  Zi)  ]  —  2, ...  ,k  . 

(53) 


«=0 


Hence 


Qiki2])  —  n  &  ~  ?«)  fl(.i  ~  z»)  ~ 

(53)  *=2 

it  k 

=  n  (**  -  z»)  ii(2i  —  z») = Qiki2}) 

(53) 

„  k  *  fc  fc 

» = n^~  *•■)  n(*_  *) = n  (*■*») 

1=0  1=2 

=  Q2*(^)  • 


Therefore, 


JVTTT-  Qui2)  1 


Using  (3.4b)  and  (3.14)  in  (3.6)  we  get 


Pu- 1  w = E  «;«/(*,) + E  w(*>)  =  ft»-i(*> 

;=2  ;=0 


(3.11) 


(3.12) 


zt)=n^ 

i=2 

i=0 

(3.13) 

In'  Jn' 

II 

H 

VI 

o 

II  VI 

•"i  IM 

(3.14a) 

-  =  **(*) 

2  <  j  <  k  . 

(3.146) 

(3.15) 


and  the  proof  is  concluded. 

Theoretically,  it  is  possible  to  write  an  algorithm  based  on  (3.5).  It  means  to  approx¬ 
imate  f(A)  by  Pm(A)  where 


Pm[A)  —  °o-f  +  aiA  •+ - h  amAm  . 


(3.16) 


This  algorithm  cannot  be  used  for  practical  problem  since  huge  roundoff  errors  result.  We 
would  like  to  stick  to  the  Newton-form  which  is  much  more  robust  from  the  roundoff  errors 
point  of  view.  For  this  purpose,  we  state  and  prove  the  following  theorem: 

Theorem  2:  Let  Pa*_i(z)  be  the  interpolating  polynomial  as  defined  in  Theorem  1,  written 
in  Newton-form 

fc-i 

Pik-ii2)  =b0  +  bl(z-  z0)  +  (*  -  zo){z  -  zx)  S,(z)P,(z) 

i=i 


(3.17) 
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Si(z)  =  bji  +  b2i+i(z  —  Zi+i) 

(3.18a) 

i 

Piiz)  =  IT (z  ”  zi){z  ~  *  =  2,. . . , k  —  1 

;=2 

(3.186) 

i2i(s)  =  1  . 

(3.18c) 

then 

k-l 

P2k-i{z)  =  bo  +  bi(z  —  2o)  +  [z  —  zo){z  -  zi)  S-i(z)Ri(z) 

«=t 

(3.19) 

where 

S,R{z)  =  6&  +  6&+1(z  -  z?+1) 

(3.20a) 

6£  =  Real  (b7i) 

(3.206) 

b2i+l  =  (&M+l)  =  &2»+l 

(3.20e) 

z,*i  =  Real  (z«+i). 

(3.20 d) 

Proof.  We  have 

^2fc+l(2)  =  ft*-l(*)  +  (z  -  Z0)(z  -  Z1)Sk(z)Rk(z). 

(3.21) 

Define 

it  -  The  set  of  all  polynomials  with  real  coefficients. 

Then,  by  theorem  1 

jP2*+l(*)>  -Pjfc-xU)  €  *•  • 


Rk{z)  G  ir  by  definition. 

Since 

S*(z)  =  [JWiW  ~  Pik- i(*)l/(*  -  *>)(*  "  zi)Rk{z) 


then 


Sk  €  ir  . 


(3.22) 


Thus 


b2k+i  —  hk+i 


(3.23) 
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fak  —  &2fc+l2*+l  —  —  ^2k+lzk+l  (3*24) 

and  the  proof  is  concluded. 

Baaed  on  (3.19),  the  vector  w 

w  =  (P2fc-iU)]v  (3.25) 

is  computed  by  the  following  algorithm; 

Algorithm  2. 

ijj  =  [bol  +  bi(A  —  zqI)\v  (3.26a) 

r  =  [A  -  zoI){A  -  zxI)v.  (3.266) 

For  *  —  1, , . . . ,  k  —  1  do: 

(3.26c) 

u  =  w  +  bj{r  +  b*+lr  (3.26 d) 

'  =  +  (*/*,  =  J„(l (+,)).  (3.26c) 


This  algorithm  requires  three  vectors.  It  is  possible  to  save  one  vector  by  using  am 
algorithm  bawed  on  (3.16).  As  mentioned  previously,  this  algorithm  has  the  disadvantage 
of  sensitivity  to  roundoff  errors.  Thus, it  can  be  used  only  with  low  degree  polynomials. 
Another  possible  way  of  saving  one  vector  amd  which  is  much  more  robust  from  the  roundoff 
errors  point  of  view  is  to  apply  (3.19)  through  calculating  the  roots  of  P2*-i( *)•  Since 
Pik-ii*)  la  a  polynomial  with  real  coefficients,  we  have  the  following  set  of  roots 

•  •  •  >  Mi » •  •  •>  M*  (3.27) 


(3.30) 
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Operating  with  (3.29)  requires  only  two  vectors.  This  approach  is  limited  mainly  by  the 
sensitivity  of  the  algorithm  which  finds  the  roots  of  P2 

4.  Rate  of  Convergence. 

We  have  the  following  definitions: 

Definition  (4.1):  Let  Tr  be  the  image  under  0  of  the  circle  |u>|  =  R{R  >  p)  and  Ir  is 
the  closed  Jordan  region  whose  boundary  is  Tr.  If  /(z)  is  single  valued  and  analytic  on 
Ir,  then  the  sequence  of  polynomials  Pm(z)  is  said  to  converge  to  f(z)  on  Ir  maximally  if 

\f(z)-Pm(z)\<C(p/R)m  zeiR  (4.1) 

where  C  depends  on  pj R  but  not  on  m  or  z. 

Definition  (4.2):  The  set  of  interpolating  points  z,  =  0(u/y)  is  said  to  be  uniformly 
distributed  on  Tp  (the  boundary  of  D)  if  wy  are  equally  distributed  on  the  circle  |u/|  =  p. 
Using  these  two  definitions,  we  can  quote  the  following  [Wals56]: 

Theorem:  Let  D  be  a  closed  Jordan  region.  Let  the  points  lie  on  the  boundary 
To  of  D.  A  necessary  and  sufficient  condition  that  the  sequence  of  polynomials  Pm{z) 
of  degree  m  found  by  interpolation  to  a  function  f(z)  analytic  on  D  in  the  points  0^ 
converges  uniformly  to  f(z )  on  D  is  that  the  points  0^  be  uniformly  distributed  on  To- 
If  this  condition  is  satisfied,  the  sequence  Pm(z)  converges  maximally  to  f(z)  on  D. 

Given  a  domain  D  and  a  function  f{z),  p  and  R  are  defined  explicitly  and  we  have 

pjR  —  the  asymptotic  rate  of  convergence. 


If  /(a)  is  an  entire  function,  (4.1)  is  satisfied  for  arbitrary  R.  Using  Theorem  (1), 
we  can  expect  that  the  algorithm  described  in  Section  3  will  converge  very  rapidly  for 
the  approximation  of  the  operator  exp  (A).  On  the  other  end,  for  the  operator  A-1,  the 
rate  of  convergence  will  depend  strongly  on  the  size  of  D  (the  parameter  p)  and  on  its 
orientation  with  regard  to  the  singular  point  z  =  0  (the  parameter  R  —  j<£(0)j).  Since  the 
set  of  interpolating  points  Zj  =  0  )  j  =  0, ...,m  —  1  depend  on  m,  given  a 
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desired  accuracy  ff,  one  has  to  decide  a  priory  on  the  degree  of  the  polynomial.  Deciding 
on  m  can  be  done  in  the  following  ways: 

1.  Getting  the  parameters  p  and  R  (analytically  or  numerically)  and  choosing  m  such 
that 

{p/R)m  S  c  .  (4.2) 

2.  Computing  the  error  numerically  on  a  set  of  check  points  on  the  boundary  for  different 
m’s,  and  choosing  m  which  will  satisfy  the  desired  accuracy. 

Using  1  or  2  gives  us  information  on  j f(z)  -  Pm(«)|00.  Substituting  it  for  ||/(^4)  - 
Pm  (A)  ||  oo  is  relatively  accurate  when  A  is  a  normal  matrix,  since  in  this  case  we  have 

ii  i(a)  -  p„(a)iu  <  i/w  -  Pmwuimr'iini 

=  I  f[z)  ~  -Pm(2)|oo 

where  T  and  T~l  are  the  unitary  matrices  which  diagnolize  A.  When  A  is  “far”  from 
normality,  m  should  be  increased  by  an  amount  which  depends  on  ||T||  •  ||T,-1||. 

5.  Applications. 

Frequently,  while  solving  a  system  of  O.D.E.’s  or  P.D.E.’s,  we  end  up  with  the  following 
set  of  differential  equations 

d 

jjUff  —  GnUn  =  Sfif(x,t)  (5.1) 

tMo) 

where  U s  and  Ss  are  vectors  of  dimension  N  and  Gs  is  a  N  x  N  matrix.  Expanding 
Ss(z,t)  as 

k 

Ss{x,t)  =  '%2aj{t)S’N{x) 

>=i 

enables  us  to  write  the  formal  solution  of  (5.1)  as 

* 

UN{t)  =  exp  {GNt)U°N  +  £  f}(Gst)S3N(x) 

y=i 


(5.2) 
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where 

//(Gy*)  =  /  exp(GNT)ctj(t  -r)d  .r  1  <  j  <  k  (5.3) 

Jo 

Us(t)  can  be  approximated  by  the  algorithm  discussed  in  Section  3  where  the  functions 
which  have  to  be  interpolated  are 


fo(z)  —  exp(tz)  (5.4  a) 

fj(z)  =  f  exp(zr)a,(t  -  r)dr  1  <  j  <  k  (5.46) 

J  o' 

In  the  case  where  (5.1)  is  originated  from  a  system  of  hyperbolic  P.D.E.’s,  the  domain  D 
is  on  the  imaginary  axis  (in  the  constant  coefficients  case)  or  close  to  it  (in  the  variable 
coefficients  case).  When  (5.1)  is  originated  from  a  set  of  parabolic  equations  the  domain 
D  is  on  the  negative  real  axis  or  close  to  it.  In  these  two  cases,  the  Faber  polynomials  are 
scaled  and  translated  Chebyshev  polynomials.  Thus,  an  efficient  algorithm,  which  makes 
use  of  the  three  terms  recurrence  relations,  can  be  implemented.  These  two  cases  are 
described  and  treated  in  details  in  [Tale86j,  [KoTa86],  [TaKo84],  [Tale85]. 

In  the  more  general  situation,  when  we  have  both  parabolic  and  hyperbolic  terms  in 
the  equation,  the  domain  D  is  more  complicated.  Let  us  look  at  the  following  simple  1  —  D 
equation 

ut  =  aw  +  6ux  4-  cuxz  [a  <  0,c  >  0)  .  (5.5) 

If  the  solution  is  periodic  in  space,  then  a  good  approximation  can  be  achieved  by  pseu- 
dospectral  Fourier  descretization.  We  get  the  following  semidiscrete  representation  of  (5.5) 

[Us)t  =  a  Us  +  6(17jv)s  +  c(Un)zx  (5.6) 

where 

N 

uN=  Y,  (5-7) 

k——N 

is  the  projected  solution.  (5.6)  can  be  written  as 


[Us)t  =  GnUn 


(5.8) 
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while  Gff  is  the  finite  domain  operator 

Cw  =  0/  +  6|.  +  e^.  (5.9) 

The  function  e*fc*  is  an  eigenfunction  of  Gs-  Thus,  if  A*  is  an  eigenvalue  of  Gff ,  then 

A*  =  a  -  c*2  +  ibk  |*|  <  N  .  (5.10) 


Since  a  <  0  and  c  >  0,  we  get  that  the  domain  D  is  the  following  parabola  in  the  complex 
plane: 

D  =  {z  —  x  +  t'yjx  =  o  -  c*2;  y  =  bk  j*|  <  N}  .  (5.11) 

In  real  applications  a,  b ,  and  e  are  not  constant,  but  have  space  dependence.  There¬ 
fore,  the  domain  D  will  vary  accordingly.  In  section  7  we  report  on  some  numerical 
experiments  treating  this  problem. 

In  a  joint  work  with  Dr.  Dan  Kosloff  and  his  colleagues,  we  investigate  the  implemen¬ 
tation  of  the  present  algorithm  to  some  real  geophysical  problems  which  can  be  represented 
as  (5.1).  The  eigenvalues  are  scattered  close  to  a  T  shape  domain  D 

D  =  {z\zeDiUD2  ;  Di  =  [— a,0)  ;  D2  =  [~ib,ib]} 


In  this  case,  we  have  an  analytic  expression  for  0(w).  (Thanks  to  Nick  Trefethen.)  The 
conformal  mapping  0(w)  which  maps  the  complement  of  the  unit  disc  on  the  complement 


of  D  is: 


(1  +  E)(w  +  — )  -4-1  —  E  ~  4 


M  =  l . 


(5.12) 


E  =Va2  +  b2/b 


Hence 


V’(w)  =  i>(v/p) 


(5.13) 


where 

(5.14) 

is  the  logarithmic  capacity  of  D.  Numerical  results  for  this  domain  are  presented  in  section 


7. 
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A  very  important  area  where  our  method  can  be  implemented  is  the  one  of  solving  general 
nonsyznmetric  systems  of  equations.  This  is  topic  of  the  next  section. 

6.  Linear  System. 

The  iterative  solution  to  a  set  of  linear  equations  which  can  be  written  in  matrix  form 

Ax  =  b  (6.1) 

is  a  well  treated  problem  in  the  numerical  analysis  literature.  While  very  efficient  algo¬ 
rithms  have  been  developed  for  the  case  when  A  is  a  symmetric  positive  definite  matrix, 
the  general  nonsymmetric  case  is  still  a  challenging  problem. 

In  this  section,  we  would  describe  an  iterative  algorithm  for  the  solution  of  (6.1)  based 
on  the  new  approach.  Since 

x  =  A~lb  (6.2) 

we  can  write  the  numerical  approximation  Xk  as 

x*  =  Pk(A)b  (6.3) 

where  Pk{z)  is  “almost  best”  approximation  to  the  function  j.  z  belongs  to  the  domain 
D  which  includes  all  the  eigenvalues  of  A. 

In  [Mant77],  (Mant78),  T.A.  Manteuffel  has  described  an  iterative  algorithm  based  on 
enclosing  the  eigenvalues  in  ellipses  in  the  complex  plane,  thus  getting  an  approximation 
based  on  scaled  and  translated  Chebyshev  polynomials  expansion.  The  algorithm  described 
in  our  paper  is  more  general  since  it  can  be  implemented  to  any  domain  in  the  complex 
plane.  Its  advantage  will  be  significant  when  the  discrepancy  between  the  best  ellipse  as 
defined  in  [Mant77]  and  the  domain  is  relatively  large. 

A  standard  approach  for  solving  (6.1)  is  known  as  Richardson  algorithm.  It  takes  the 
following  form: 

xk+l  =xk_  ak(Axk  _  j)  (6 .4) 

If 

Ek  =  x-xk  (6.5) 
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we  get 

£*  =  |P;(A)|B°  (6.6) 

where 

fc-i 

P»'(A)=  (6.7) 

i=o 

and  ay  are  optimal  in  the  sense  that 

majc j Pk (*)|  <  nwx \Pk{z) |  V  Pk  €  irk  (6.8) 

*€£>  *6  D  f»k(0)=l 

(7T*  is  the  space  of  all  polynomials  of  degree  k).  We  would  like  to  show  that  am  algorithm 
which  results  from  implementating  our  approach  is  equivalent  to  algorithm  based  on  (6.4)- 
(6.8).  For  this  purpose  we  have 

Theorem  3.  Let  Tk(z)  be  a  polynomial  of  degree  k  defined  onC  such  that  Tfc(O)  =  1  and 
let  Qk~i(z)  be  the  inter polant  of  the  function  j  at  the  roots  of  Tk  then. 


1  -  a  g*_x(z)  =  Tk(z)  .  (6.9) 

Proof.  We  have 

Qk-i{zi)  =  -  *  =  1 . k  (6.10) 

Z\ 

where  Zi  is  a  root  of  Tk.  Therefore,  the  polynomial  Rk(z) 

Rk(z)  =  Tk(z)  -  (1  -  z  Qk-i[z))  (6.11) 


vanishes  at  k  +  1  points:  0,*i, . . . ,  zk  and  thus  it  is  identical  zero.  Hence, 

Tk(z)  =  1  -  z  Qfc_x(z)  (6.12) 

and  the  proof  is  concluded. 

If  x°  is  the  initial  guess,  we  have 

Ax°  =  b°  (6.13) 

A( x  -x°)  =  b-b°  (6.14) 
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x  =  A-l(b-b°)+xQ  . 

(6.15) 

Approximating  A-1  by  the  interpolant  <3*_i(A)  results  in 

*k  =  Qk-i(A)(b-b°)+x° 

(6.16) 

using  (6.5)  and  (6.14)  we  get 

Ek  =  [I  -  AQk-X{A)\E°  . 

(6.17) 

Since  (6.9),  the  equivalence  is  established. 

We  have  shown  that  using  the  algorithm  based  on  interpolating  the  function  j  at 
will  reproduce  identical  results  to  Richardson  iterations 

xJ+1  =  x*  -  aj{Ax*  -  6) 

(6.18) 

with 

<*j  -  —  j  =  o, . . . ,  k  -  1  . 

zi 

(6.19) 

Writing  (6.18),  (6.19)  as  an  algorithm  with  real  arithmetic  takes  the  following  form: 
Algorithm  3  (Richardson). 

(The  parameters  are  ay  =  l./zy  and  zy  are  defined  in  Appendix  A.) 

x1  =  x°  -  ao(Ai°  -  b) 

(6.20a) 

x2  =  x1  ~  a^Ax1  -  b) 

(6.206) 

for  *  =  1,  ...  , 

R*  =  Ax2i~l  -  b 

(6.21a) 

x2‘+l  =  x2‘“l  -  [2a?  -  |aij3A)i2*  (of  =  Re(a<)) 

(6.216) 

Preconditioning. 

Usually,  solving  a  linear  system  of  equations  (6.1)  is  composed  of  two  stages: 
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1.  Modifying  the  original  system  to  an  equivalent  one 

Ax  =  b 


(6.22) 


such  that  (6.22)  will  converge  faster  than  (6.1). 

2.  Solving  (6.22). 

In  many  cases,  we  have  a  family  of  matrices  A  which  depend  on  a  parameter  u  such 

that 

Aux  =  L  (6.23) 

and  we  would  like  to  choose  the  optimal  u.  Based  on  Section  4,  it  is  obvious  that  the 
significant  factor  which  determines  the  convergence  is 

\  =  p/R.  (6.24) 

Optimization  is  achieved  when 

|AK*)|  <  |A(w)|  Vu;  .  (6.25) 

Therefore,  we  conclude  that  for  matrices  whose  eigenvalues  are  scattered  in  the  complex 
plane,  one  should  consider  the  factor  p/R  rather  than  the  standard  condition  number 
cond(A)  =  ||A||  ||A-1||.  It  is  possible  to  have  two  matrices  A  and  B  such  that 

cond  (A)  <  cond  (J3)  (6.26) 

and  on  the  other  hand 

(p/fip)  >  WR\{B)  .  (6.27) 

For  example:  Let  A  be  a  normal  matrice  whose  spectrum  r(A)  is 

r(A)  ={z\l<  |*|  <  2;  Re  z  >  0}  . 

Let  B  be  a  normal  matrice  whose  spectrum  is 


r{B)  =  {*|  |*  -  3|  <  2}  . 
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We  have 

cond(A)  =  2 
con d(B)  =  5  . 

On  the  other  hand 

\plR\(A)  =  0.8547 
[p/R\{B)  =  0.6666*  ••  . 


(6.28) 


(6.29) 


Since  (6.29),  Richardson  algorithm  for  B  will  converge  much  faster  than  for  A,  even  so 
coad(B)  >  coad(A).  Based  on  this  discussion,  we  would  like  to  show  also  that  the  standard 
S.O.R.  procedure  can  be  considered  as  an  incomplete  optimization  process.  Solving  (6.1) 
by  S.O.R.,  we  use  the  following  iterative  procedure 


x*+1  =  Twxk  +  bw 


(6.30) 


where 

Tu  =  MZ'NU 
K  =  AC1* 

Afw  =s  D  +  ujli 

Nu  =  (1  -  u)D  -  uU  . 

U,  Dy  L  are  the  lower,  diagonal,  and  upper  parts  of  the  matrice  A  respectively.  The 
optimal  wior  is  chosen  such  that 


|r(rw.J|  =  minimal  . 


(6.32) 


(r(A)  is  the  spectral  radius). 

If  x  is  a  solution  of  (6.30),  then 

Awx  =  bu  (6.33) 


where 

Au  =  I  -  T„  .  (6.34) 

Thus,  it  is  evident  that  rather  than  using  an  optimization  procedure  based  on  (6.32) 
one  should  use  the  more  general  one,  based  on  (6.25) .  A  well-known  example  treated  in 
the  literature  for  demonstrating  the  features  of  S.O.R.  is  the  problem  of  solving  Laplace 
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equation  in  a  rectangle.  In  (Youn7l,  p.  205],  it  is  shown  that  for  a  certain  discretization 
of  the  rectangle,  the  optimal  u  is 

Uior  =  1-25  • 

In  this  case,  all  the  eigenvalues  A*  of  Tu  satisfy 

| A,- 1  =  0.25  . 

Hence,  the  rate  of  convergence  is  0.25.  Optimizing  with  respect  to  (6.25)  can  result  with 
a  different  solution.  We  do  not  intend  to  carry  out  this  optimization  process,  but  to  point 
out  a  different  possible  parameter  u>.  For  this  problem,  we  can  choose  u 

U  =  1  +  e  0  <  e  <  0.25 

such  that  the  domain  D  which  includes  ail  the  eigenvalues  of  A&  will  be 

D  —  Bi  U  B2 


where 


and 


B  i  =  <A0} 

B2  =  {z|  \z-l\<  e} 
0  <  Aq  <  1  -  s  . 


Since  the  complement  of  D  is  not  a  simply  connected  domain,  it  is  not  included  in  the 
theory  discussed  in  this  paper.  Implementing  our  approach  to  these  types  of  domains  will 
be  carried  out  in  a  future  paper.  For  the  time  being,  let  us  mention  that  the  basic  results 
extend.  Thus,  since  B\  is  composed  of  only  one  point  and  B2  is  a  circle  of  radius  s  around 
1,  it  can  be  shown  that  the  rate  of  convergence  p/R  is 


p/R(A&)  =  e  <  0.25 


and  Ao  will  enter  into  the  constant  C  (4.1).  The  set  of  interpolating  points  is  the  union  of 
sets  of  both  domains.  Since  in  B i  we  have  only  one  point,  we  will  get  the  following  points 


^0,  Ml*  *  •  •  *  Mn 
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while  (ij  can  be  chosen  as  equally  distributed  on  the  circle  \z  -  1\  =  e 

Hj  =  s  exp(2irtj/n)  j  =  1, . . . ,  n 

or  as  the  n  zeros  of  Faber  polynomial  of  degree  m  which  is  ( z  —  l)n;  thus 

Mj  =  1  j  =  l,...,n. 

The  efficiency  of  using  an  algorithm  based  on  w  rather  then  u/Jor  depends  on  the  accuracy 
needed,  since  the  constant  C  has  been  increased  by  a  factor  of  l/Ao. 

7.  Numerical  Results. 

I.  Time  dependent  problem- Parabolic  type. 

In  this  subsection  we  present  numerical  results  for  the  following  problem 

u  t  —  Gu  =  S{x,t)  (7-la) 

u(x,0)  =  0  . 

Where 

G  =  a(x)  +  6(i)  —  +  c(x)  (7.16) 

S(x,t)  =  sin(fcx)  -  t  {(-fc2a(x)  +  c(x)]  sinfcx  +  kb(x)  cosfcx}  .  (7-lc) 

The  exact  solution  of  (7.1)  is 


u(x,  t)  =  t  si nkx  . 

In  order  to  solve  (7.1)  numerically,  we  first  approximate 
difference  operator  Gn .  Assuming  periodicity  we  have 


{GNu)j  =  o(x>) 


u j+i  -  2 u,  +  Uj-i 

Ax2 


+  bixi) 


*»;  + 1  ~  U3  - 

2Ax 


where  N  i is  the  number  of  grid  points  and 


(7.2) 

the  space  operator  G,  by  the  finite 

—  ■+■  c(xy)uy  J=0 . IV -1 

(7.3) 


Ax  =  2ir/N 


(7.4) 
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x j  =  j  •  Ax 


j  =  0,...,JV  . 


Hence  Gpt  ia  a  N  x  N  matrix.  The  semidescrete  representation  of  (7.1)  is 


3  =  0 . iV 


(uw)t  -  Gsun  =  Sff  +  tSft 

K)i  =  o  y  =  o, . . . ,  jv 

(5^)j  =  sin  (Axy)  j  =  0, ...  ,N 

(Stf)j  =  [k?a(xj)  -  e(xy)j  sin  kxj  -  kb(x})  coe  Axy  j  =  0, 

The  formal  solution  of  (7.6)  is 

uiv(t)  =  fi(G^t)S^  +  /2  [Gfft)Sjf 


(7-6 a) 
(7.66) 
(7.6c) 
.,AT.  (7.6d) 


where 


WUCIC 

/i(GatO  =  f  exv(GNr)dr  =  (exp(GVt)  -  I\/Gn  (7.8 a) 

Jo 

f2{Gst)  =  f  exv(GNr)(t  -  r)dr  =  (exp(GiVf)  -  Gy*  -  /]/G*.  (7.86) 

Jo 

In  order  to  implement  our  algorithm  we  have  to  get  an  approximation  of  the  domain  D 
which  includes  the  eigenvalues  of  Gs-  One  way  to  get  it  is  by  doing  a  Fourier  analysis  of 
the  constant  coefficients  operator. 


(Gjvti)y  =  a 


_Uy+X  -2uy  +  Uy— i  tty+i  ~  uy-i  _ 


(7.9a) 


where 


:  max  |a(z)|;  6=  max  |6(x)l; 

0<*<a»'  V  71  0<x<2ir  v  n 


c  =  o3SJc(l)l!  or  c  =  ■ 

(7.96) 


Let  be  an  eigenvector  of  G s',  then 


(Vfc)y  =  « 


_  iky  Ax 


(7.10) 


A*  =  ^j-(coe(*Ax)  -  1)  -  *^sin(fcAx)  -  c 


(7.11) 
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is  an  eigenvalue  of  Gs. 

The  tables  in  this  subsection  present  numerical  results  for  the  fully  descrete  solution 
of  (7.1)  where 


k  =  3 


a(x)  =  l./(2  +  coax) 

(7.12) 

b(x)  =  1-/(2  +  sin(x)) 

(7.1?) 

c(x)  =  — 20./(2  +  cosx). 

(7.14) 

From  (7.11)  we  get 

A<Re\k<B 

(7.15) 

\Im*k\<C 

(7.16) 

while 

A  = 

--1.-20;  B  =  -20/3  ;  C=±. 

(7.17) 

Since 

\B-A\^>C 

(7.18) 

taking  D  as 

D  =  {z(A  <  x  <  B} 

(7.19) 

will  be  a  relatively  good  approximation  to  the  domain  D .  We  have 

p(D)  =  {B-A)/ 4. 

(7.20) 

Hence,  in  order  to  implement  the  new  algorithm  we  have  to  interpolate  the  functions 

fi{pzt)  =  [exp [pit)  -  1  }/pz 
h{pzt)  =  [exp {pit)  -pH-  1  }/(pz)2 

at  points  on  jD. 

In  this  case  we  can  take  the  M  interpolating  points  as  the  extreme  of  the  scaled  and 
translated  Chebychev  polynomial  of  degree  M  [Rivl74).  Thus 

5j  =  ^[(B  —  A)ij  +  B  +  A]  j  = 
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where 


ij  =  coe 


O'  -  *)*' 

M-  1 


In  order  to  avoid  roundoff  errors  one  should  arrange  the  interpolating  points  as  described 
in  Appendix  A.  Hence 


where 


zi  ~  ~  ^)si  + 


j  = 


(7.21) 


Xi  =  1  ;  x2  =  -1 


(7.22a) 


—  -R«(w2;-3)  j  =3 ,. . .  ,M  (7.226) 

and  way-3  are  as  in  Figure  (A.l)  for  the  case  M  =  7. 

Since  zy  are  real  we  can  use  algorithm  1  (3.2).  The  next  two  tables  present  the 
numerical  results  using  algorithm  1  with  zy  defined  by  (7.21)-(7.22). 

Table  1. 

Mesh  refinement  chart  -  problem  (7.1) 

Using  algorithm  1.  t—  1 


N 

M 

L2In 

32 

24 

.1108-01 

64 

52 

.2592-02 

128 

112 

.7033-03 

N  -  Number  of  grid  points. 

M  -  Number  of  matrix-vector  multiplications  (Each  evolution  operator  is  approximated 
by  a  polynomial  of  degree  M/2) 

La/„  -  L2  Error  at  t  =  1. 

For  sake  of  comparison  we  present  in  Table  2  a  similar  chart  while  using  a  standard 
second  order  in  time  scheme 

A  #2 

un+l  =un  +  AtGNun  + —G2Nun  (t  =  n  •  At). 


(7.23) 


Table  2. 

Mesh  refinement  chart  -  problem  (7.1) 
Using  (7.  SS)  t  =  1 


N 

M 

L3Ta 

32 

102 

.1108-01 

64 

408 

.272^-02 

128 

1632 

.6828-03 

(In  this  case  M/2  is  the  number  of  time  steps.) 

Comparing  Tables  1  and  2  we  observe  that  while  in  the  standard  second  order  scheme 
(7.23)  M  is  proportional  to  N2  (which  results  from  the  fact  that  At  =  0(Ax2))  for  the  new 
algorithm,  M  is  “almost”  proportional  to  N.  This  phenomenon  is  explained  and  proved 
in  [Tale  85).  For  t  =  20  we  have  the  following  results; 

Table  3. 

Mesh  refinement  chart  -  problem  (7.1) 

Using  algorithm  1  t  =  20 


N 

M 

L2In 

32 

24 

.1322-01 

64 

52 

.3407-02 

128 

112 

.1180-02 

Observing  Tables  1  and  3  we  notice  that  M  does  not  depend  on  t.  It  can  be  explained 
as  follows:  For  large  t,  ||  exp(G//t)||  is  much  smaller  than  the  error  which  results  from  the 
space  descretization.  Thus,  since  (7.8)  we  have 

£  -1  (Gff 

h{GNt)~-{GNt  +  I)lG2N 

which  means  that  for  large  t,  approximating  fi(Gst)  is  equivalent  to  inverting  Gjv.and 
approximating  f{Gst)  is  equivalent  to  inverting  G2N. 
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We  did  not  implement  the  second  order  algorithm  (7.27)  for  t  =  20,  but  it  is  straightforward 
(by  stability  considerations)  that  for  N  =  32, 64, 128  we  have  M  =  2040, 8160, 32640 
respectively.  Hence  for  t  =  20  and  N  =  128  the  new  algorithm  is  more  efficient  than  the 
second  order  scheme  by  the  impressive  factor  of  32640/112  =  290. 


II.  Viscoelastic  Wave  Propogation. 

We  have  worked  on  this  model  with  the  Geophysical  group  in  Tel-  Aviv  University 
headed  by  Dr.  Dan  Kosloff.  A  detailed  report  on  the  results  for  a  general  2-D  problems 
will  be  published  elsewhere.  In  this  subsection  we  describe  the  implementation  of  the 
algorithm  for  the  simple  1-D  model 


where 


Ut  =  GU  0  <X<L 


(  ° 

1 

0  \ 

(Ul\ 

G  = 

m„V 

0 

V 

;  u  = 

u2 

\  a 

0 

—  l/r  ) 

Vu3  / 

V  = 


dx 2 


mu 


1  -  (1  - 


c  =  2000  (the  speed  of  sound) 


a  =  c2(l  -  r/rj)/ri 


(7.24) 


(7.25) 

(7.26a) 

(7.266) 

(7.26c) 

(7.26d) 


xi \  is  the  pressure, u2  is  the  pressure  time  derivative, U3  is  a  memory  variable  and  r,  >7  are 
parameters  of  the  problem.  The  initial  data  are 


ui(z,0)  =  exp 


(7.27a) 


u2(z,0)  =  u3(z,0)  =  0. 


(7.276) 
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The  space  descretisation  is  done  by  Fourier  method  (GoOr81).  Thus  the  semidescrete 
representation  is 

(^w)t  =  GffUu  (7.28a) 

=  P«[«p(-i(i  -  ii)Ji .  (?284) 

uUx.o)  =  crJ(i,o)  =  o . 

Ps  is  the  pseudospectrai  projection  operator: 

s 

Pn/  =  5Z  ak  exp(»(23r*/L)*)  (7.29) 

fc=-jv 


= 


1 

2iVC* 


2JV-1 

53  /(*y)  exp(-»(2jrfc/L)x,] 


,=o 


C*  =  1  for  |fc|  iV 
C/v  —  C_iv  =  2 


(7.30) 


and 


=  j  Ax 

j -0,...,2N -1 

Ax  =  L/2N 


Gff  =  PffGPs  • 


[Gs  is  a  6 N  x  6 N  matrix.) 

In  our  experiments  we  took 

N  =  64  . 


Thus 


Ax  =  20  . 


L  =  2N  x  Ax  =  2560  . 


(7.31) 

(7.32) 

(7.33) 

(7.34) 

(7.35) 


Since  N  =  64,  G#  is  a  384  x  384  matrix. 

Approximating  the  domain  D; which  includes  the  eigenvalues  of  Gjy.is  done  by  making 
use  of  the  following  idea.  Let  us  assume  that  the  number  of  points  K  which  are  needed  to 
resolve  the  coefficients  of  the  operator  G  are  relatively  small.  In  this  case,  decreasing  the 
space  domain  by  a  factor  of  K/N  and  using  the  same  Ax  gives  us  the  matrix  G*.  Since 


K  is  small,  one  can  use  a  library  routine  to  compute  the  eigenvalues  of  Gk-  The  domain 
Dk, which  includes  this  set  of  eigenvalues, is  relatively  a  good  approximation  of  D.  In  the 
present  case  Gff  is  a  constant  coefficient  matrix, and  the  domain  D+  achieved  by  computing 
the  eigenvalues  of  G*  (a  12  x  12  matrix)  is  exactly  the  same  as  Du%.  The  domain  D4  is 

£>4  =  {z\z  €  [-o,0]  or  z  €  [ — *6,  *6}}  (7.36) 

where  a  and  6  depend  on  r,  rj.  For  this  domain  we  have  an  analytic  expression  of  the 
conformal  mapping  which  maps  the  complement  of  the  unit  disc  on  the  complement  of  D\ 
(5.12).  The  logarithmic  capacity  is  given  by  (5.14). 

We  ran  two  sets  of  numerical  experiments.  In  the  first  one  we  took 

r  =  0.1600890  x  10“3 


and  we  have  gotten 


V  =  0.1582262  x  10“3 


a  =  6320  ;  b  =  317  . 


Using  algorithm  2  (3.26)  for  computing  the  solution  at  time  level  t  =  0.1  results  in 


N 

M 

£2/n 

128 

130 

.1588-02 

128 

150 

.2904-05 

(Since  we  do  not  have  am  analytic  expression  for  the  exact  solution,  we  computed 
£>3£»  by  comparing  the  numerical  solution  to  another  numerical  solution  achieved  by  using 
algorithm  2  with  M  =  3 00.)  Similarly,  using  the  second  order  in  time  algorithm  (7.23) 
results  in 


N 

M 

L7Ta 

128 

632 

.1902-1 

(For  M  <  632  the  scheme  is  unstable.) 
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In  the  next  experiment  we  took 

r  =  0.1600890  x  10-4 


T)  =  0.1582262  x  10~4 


and  we  have  gotten 

a  =  63201 


b  =  317  . 


The  results  at  t  =  .1  were 


N 

M 

L2ln 

128 

400 

.2200-02 

128 

450 

.1096-04 

While  using  (7.23)  results  in 


N 

M 

L2Ta 

128 

6320 

.1892-3 

(For  M  <  6320  the  scheme  is  unstable.) 


m.  Linear  Systems. 

For  the  numerical  experiments  reported  in  this  subsection,  we  have  used  a  test  matrix 
AHxlf  which  is  block  diagonal.  Each  block  is  of  the  following  shape 


Hence,  the  eigenvalues  are 


(7.37) 


\j  —  i  i |hyCyj  j  —  1, . . . ,  N/2  . 


(7.38) 


i 
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Two  kinds  of  experiments  have  been  carried  out: 
bj  =  Cj  =>  the  matrix  is  normal. 
bj  ^  Cj  =>  the  matrix  is  not  normal. 

We  solve 

Av  =  w  (.A  is  a  1000  x  1000  matrix) 

where 

V  =  [t/1,  .  .  .  ,  Uiooo]r 

and 

vj  =  {-!)’  j  =  1, . . . ,  1000 

and  the  initial  guess  v°  is  such  that 


v 


o 

i 


=  0 


J  =  l, 


,1000  . 


Three  methods  have  been  tested: 

(1)  In  -  the  method  described  in  this  paper 

(2)  Cb  -  Chebyshev  approach  (Mant78) 

(3)  Mr  -  Minimum  residual. 


The  minimum  residual  algorithm  is  defined  as  follows: 
Given  initial  guess  z°  for  k  =  0, 1, . . . ,  until  satisfied  do 

r*  =  Axk  -  b 

o*  =  (r*,Ar‘)/(Ar*,Ar‘) 


(7.39) 


(7.40) 


(7.41) 
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1.  First  case  -  Cross  shaped  domain. 

In  [SmSa85)f  the  authors  treat  an  example  given  by  Hageman  [HaYo81j  which  orig¬ 
inates  from  the  solution  of  the  neutron  diffusion  equation,  where  the  eigenvalues  of  the 
Jacobi  iteration  matrix  may  be  shown  to  lie  on  a  the  following  Cross-shaped  domain  D: 


Since  (7.42)  we  get  that 

P  -  1/V2  .  (7.43) 

Thus  the  mapping  function  from  the  exterior  of  disc  of  radius  p  to  the  exterior  of  D  is 

z  =  rl)[ w)  =  b  +  y  w2  +  .  (7.44) 

Since 

S  =  N'_l(0)|  (7.45) 


we  have 

„  \b2-h(b*-l)l^]U2 

“  L - 2 - .  • 

Using  (7.43)  and  (7.46)  we  get  that  the  asymptotic  rate  of  convergence  is 

e/R  =  [i1 . 


(7.46) 

(7.47) 


According  to  the  theory  (4.1) 


1/(4)  -  />„(*)!  <  CMR)M  ■ 


(7.48) 
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thus  we  can  predict  that  after  M  iterations,  the  accuracy  (in  the  normal  case)  will  be 

L3Error(/n)  £  {p/R)M  .  (7.49) 

An  improved  prediction,  which  includes  the  constant  C  is: 

M 

L2Error(Jn)  3  C(p/i2)"  3  IJ(1  -  -)  (7.50) 

«=i  * 

where  z  is  any  point  at  the  domain  D  and  a,  are  the  interpolating  points. 

(In  the  next  two  tables  we  have  chosen  z  =  (6,0).) 

The  numerical  results  are 


Table  7. 

Solution  of  ( 7.89 )  with  b  =  1.004 


■a 

(<>/*)" 

CWR)% 

L2  Error  (Jn) 

£,3Error  (C&) 

LaErtor  (Afr) 

Normal 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

102 

1.575-3 

3.16-3 

3.87-3 

1.02+0 

1.85-1 

9.32+0 

9.67-2 

6.55+0 

201 

2.8-6 

5.63-6 

6.9-6 

4.41-3 

1.04-1 

10.31+0 

5.37-2 

7.34+0 

402 

9.0*12 

1.8-11 

2.2-11 

2.29-8 

3.9-2 

7.65+0 

1.9-2 

7.13+0 
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Table  8. 

Solution  of  (7.S9)  with  b  =  1.1 


M 

(/>/«)“ 

C{p/R)% 

Li  Error  (7n) 

I^Error  (C(,) 

Li  Error  (Mr) 

Normal 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

10 

4.13-2 

8.99-2 

1.10-1 

9.69-1 

1.9-1 

9.9-1 

1.57-1 

8.0-1 

22 

9.03-4 

1.8-3 

2.22-3 

6.37-2 

5.01-2 

5.56-1 

4.07-2 

4.29-1 

42 

1.5-6 

3.1-6 

3.8-6 

2.82-4 

6.34-3 

1.33-1 

9.98-2 

82 

4.5-12 

9.0-12 

1.1-11 

2.15-9 

1.2-4 

4.81-3 

3.27-3  . 

2.  Second  Case  -  Boomerang  shape  domain. 

Another  example  reported  in  [SmSa85]  is  Van  der  Vorst’s  example.  Let  M  be  the 
matrix  arising  from  the  descretization  of  the  P.D.E.  operator  urI  +  «y»  +  0\uz  +  and 
let  K  be  the  incomplete  Cholesfcey  factorization;  then  the  eigenvalues  of  the  preconditioned 
matrix  K~XM,  are  sometimes  observed  to  form  the  following  boomerang-shaped  profile 


Since,  in  this  case  we  do  not  have  an  explicit  expression  for  the  conformal  mapping, 
we  have  used  numerical  algorithm  (Tref8l)  for  computing  the  interpolating  points  z,.  The 
predicted  accuracy  is 

C(o/K)S  a  II<1 " 


(7.51) 
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while 

*  =  (6,0).  (7.52) 

Table  9  summarizes  the  results  of  this  case. 


Table  9. 


M 

C(p/R)“ 

£2  Error  (Jn) 

£2  Error  ( Cb ) 

I-2Error  (Mr) 

Normal 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

20 

1.23-4 

5.4-4 

8.8-4 

isa 

1.33-1 

2.42-1 

2.35+0 

40 

7.58-8 

2.36-7 

5.14-7 

MM 

5.8-2 

1.13-1 

2.35+0 

60 

3.78-11 

1.03-10 

3.05-10 

1.71-2 

- J 

4.24-2 

5.59-2  1 

2.35+0 

In  the  next  experiment  we  have  shifted  the  domain  to  the  left. 


According  to  the  theory,  the  two  methods:  C\>  and  Mr  would  not  converge  in  this  case. 
Table  10  verified  this  fact. 
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Table  10. 


M 

C(tf/R)“ 

L<x Error  (/n) 

Lj  Error  (C&) 

LjErr 

or  (Mr) 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

20 

1.56-3 

6.65-3 

1.08-2 

1.15+0 

1.31+0 

7.09-1 

3.62+0 

40 

1.22-5 

3.8-5 

8.25-5 

4.93+0 

6.67+0 

6.76-1 

3.62+0 

60 

7.91-8 

2.16-7 

6.38-7 

23.0+0 

56.9+0 

6.58-1 

3.62+0 

3.  Third  Case  -  disjoint  intervals. 

As  mentioned  previously  (Section  2),  since  the  complement  of  D  is  not  simply  con¬ 
nected  anymore  one  should  use  a  slightly  different  theory.  In  a  future  paper  we  will  carry 
out  a  detailed  analysis. 

In  this  subsection  we  report  on  a  few  experiments  where  D  is 


The  interpolating  points  are 


D  ~  I\  U  /3 

h  —  (<*1,0:3)  ;  /a  =  [<*3,04] 

an  <  a%  <  03  <  at  . 

zl  tz2’  ■  •  • »  ZN\ »  z\>  z2*  *  •  ■  >  "Va 


(7.53) 

(7.54) 

(7.55) 

(7.56) 

*  =  (7.57) 

*  =  (7.58) 


where 


ir(i  -  1)  , 

(a3  -  at)  cos  _-y  +  <*2  +  «i 


,  1  f,  ^  "■(*'  ~  1) 

4  =  2  [(a< "  aj)  c“ 


+  as  +  a* 
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The  error  polynomial  Pn  (z)  is 


Pn(z)  =  PNl  (*)  •  Pjv,  (*)  (AT  =  Ni  +  AT,) 

(7.59) 

while 

Ni 

pnAz)  =  II(1_  */z«) 

i=l 

(7.60) 

pNi  (2) = uu  ~  z/z*)  • 

i=i 

(7.61) 

We  have 

1-PiVi  (•*)!  ~  (Pi/pi)*1  z  €  I\  . 

(7.62) 

Since  the  conformal  mapping  from  the  exterior  of  a  unit  disc  to  the  exterior  of  an  interval 

[a,b]  is 

if  l(^)  =  \  ^(6  -a){w  +  i)  +  64-0  , 

(7.63) 

we  get 

II 

1 

(7.64) 

R  =  (a  +  b  +  Vab)/ 4  . 

(7.65) 

Thus 

IJV.01  ~  (,./«.)"■  =  (^+^) 

(7.66) 

Ps3  (*)  satisfies 

|JV,(*)|  <1  z  €  I\  . 

(7.67) 

Therefore  if 

Vv^+V^T/ 

(7.68) 

then 

I^V(*)I  <c  z  €  I\  . 

(7.69) 
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Now,  for  z  €  /j  we  have 


max  l-Py^a)!  ~  ka^1 


(7.70) 


|r».(«)l  ~  Wto)"' 

(7.71) 

where 

.  /„  _  Vs*- 

V37-VS?' 

(7.72) 

Therefore,  if 

ka4l(p2/Ri)Na  <  e 

(7.73) 

then 

|P//(s)|  <  c  z  €  h  . 

(7.74) 

Using  (7.66),  (7.68), and  (7.73), we  get  that  N2  has  to  satisfy 

*r  .  log(/»i/^i“-«)  N  .  log(l/*r) 

2  >  log (P2/R2)  1  log (P2/R2)  ■ 

(7.75) 

In  Table  (II)  we  report  on  experiments  where  aj,  aj  are  fixed,  and  we  increase  03,  04 
such  that  a*  -  03  is  constant.  According  to  (7.75),  it  is  easily  verified  that  in  this  case 
N2  — ♦  N\.  The  predicted  error  is 


Er  =  max  \Pn{z)\  =  JV(a«)  =  II ^  ~  a*lz<)  II^1  “  a*/*.?)  • 

»=i  »= 1 


(7.76) 
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Table  11. 


ai  =  1  aj  =  3  a*  =  a3  +  2 


<*s 

Si 

S2 

Er 

L2  Error  (/n) 

L2  Error  (C&) 

L2  Error  (Mr) 

Normal 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

20 

12 

12 

1.9-7 

H 

1.7-5 

6.6-5 

3.5-3 

28-3 

1.0-2 

40 

12 

12 

3.2-7 

1.9-7 

5.4-4 

9.7-4 

8.7-2 

2.2-3 

80 

12 

12 

4.4-7 

2.5-7 

1.4-4 

75-3 

1.0+0 

5.4-2 

60 

12 

12 

5.1-7 

2.9-7 

3.2-4 

4.3-2 

5.9+0 

2.8-2 

1.6-1 

320 

12 

12 

5.5-7 

ESI 

6.9-4 

7.4-2 

28.8+0 

1.2-1 

1.6-1 

In  the  next  set  of  experiments, we  also  increase  the  distance  between  a3  and  a4.  The 
results  are  reported  in  Table  12. 


Table  12. 

ai  =  1  oij  —  3  a<  =  3a3 


«  3 

Si 

S2 

Er 

L3  Error  (Jn) 

L2  Error  ( Cb ) 

L2  Error  (A/r) 

Normal 

N.  Normal 

Normal 

N.  Normal 

Normal 

N.  Normal 

10 

32 

1.6-7 

na 

1.3-5 

7.7-7 

3.6-5 

3.8-4 

2.5-4 

20 

10 

38 

8.4-7 

2.9-7 

9.7-6 

2.0-5 

1.1-3 

4.2-4 

1.5-3 

40 

10 

44 

2.2-6 

6.0-7 

9.8-6 

3.6-4 

2.2-2 

4.9-3 

1.3-2 

80 

10 

50 

3.9-6 

9.6-7 

1.5-5 

4.2-3 

2.8-1 

1.9-2 

2.5-2 

In  the  last  set  of  experiments,  we  have  negatives  eigenvalues  as  well.  In  this  case,  one 
cannot  use  the  two  methods:  C Mr. 


1 
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Table  IS. 


oti  =  2c*2  a$  =  -a?  =  —ai 


<*1 

Ni 

n2 

Er 

Li  Error  (J») 

Normal 

N.  Normal 

2 

16 

16 

9.7-7 

mfmm 

2.5-5 

4 

16 

16 

2.5-5 

8 

16 

16 

fiSSI 

II 

2.5-5 

8.  Conclusions 

It  has  been  shown  that  having  an  a  priori  knowledge  on  the  distribution  of  the  eigen¬ 
values  of  a  matrix  A,  it  is  possible  to  construct  an  efficient  algorithm  for  approximating 
f{A).  The  more  accurate  we  locate  the  domain  of  the  e.v.,  the  more  efficient  is  the  al¬ 
gorithm.  Two  methods  addressing  this  problem  were  described  in  Section  7.  It  seems 
to  us  that  once  this  problem  of  finding  the  domain  of  the  e.v.  is  solved  satisfactorily, 
the  algorithms  described  in  the  paper  can  be  used  as  a  numerical  tool  for  many  practical 
problems. 
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Appendix  A. 

The  set  of  interpolating  points  satisfy 

Zj  =  j  =  1 . N  (A.  1) 

while  Wj  are  roots  of  the  equation 

w*  =  1  .  ( A.2 ) 

The  problem  is:  how  to  arrange  w}  such  that  the  roundoff  errors  will  be  minimum. 

In  a  future  paper  we  will  carry  out  a  detailed  analysis  addressing  this  problem.  Ac* 
cording  to  the  solution  described  there^we  have 

wl,w2t...,ws  (A. 3) 

equally  distributed  on  the  unit  circle,  where 

twi  =  1  ;  u>2  =  -1 

and  each  “new”  point  Wj  is  chosen  such  that  the  partial  set 

U>l,t£>2,...,U>;  (A.4) 

will  be  as  equally  distributed  as  possible.  For  example,  when  N  =  12  we  have 


In  the  case  of  algorithm  (3.27),  (A.3)  has  to  be  modified  as  follows 


(A.5) 


ttfj,  Wj,  W$ ,  -  ,  WN 


t 
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(A.6) 

(A.7) 


Wi,W2,...,W2j-x,W2j  (A.  8) 

will  be  as  equally  distributed  as  possible.  Thus,  for  the  case  N  =  12  we  have 


Figure  (A. 2). 
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